Do Athermal Amorphous Solids Exist? 
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We study the elastic theory of amorphous sohds made of particles with finite range interactions 
in the thermodynamic limit. For the elastic theory to exist one requires all the elastic coefficients, 
linear and nonlinear, to attain a finite thermodynamic limit. We show that for such systems the 
existence of non-afhne mechanical responses results in anomalous fluctuations of all the nonlinear 
coefficients of the elastic theory. While the shear modulus exists, the first nonlinear coefficient 
B2 has anomalous fluctuations and the second nonlinear coefficient S3 and all the higher order 
coefficients (which are non-zero by symmetry) diverge in the thermodynamic limit. These results 
put a question mark on the existence of elasticity (or solidity) of amorphous solids at finite strains, 
even at zero temperature. We discuss the physical meaning of these results and propose that in 
these systems elasticity can never be decoupled from plasticity: the nonlinear response must be very 
substantially plastic. 
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I. INTRODUCTION 

Amorphous solids are ubiquitous in nature and in tech- 
nology, spanning well known materials from obsidian to 
metallic glasses. They are typically obtained by fast cool- 
ing a liquid such that it glides over the melting point 
where normally it can form a crystalline solid via a stan- 
dard first-order phase transition. Rather, the fluid be- 
comes super-cooled and experiences that so-called glass 
transition where its relaxation time becomes ever longer 
upon decreasing the temperature [ij. When T — we 
have typically a frozen amorphous structure which ap- 
pears to react elastically to small strains. In this paper 
we consider only amorphous solids that are made of TV 
point particles in d dimensions, in a volume V , which in- 
teract with each other via finite range interactions. The 
subject of our investigation is how athermal amorphous 
solids react to external strains. 

The potential energy in our strained amorphous solids 
can be written as J7({ri(7)}, 7), where {ri}^^^ are the 
positions of the particles and 7 is the applied strain. In 
this paper we shall treat the strain as a scalar, all equa- 
tions can be written in full tensorial form if required (cf. 
0), but for our purposes a single scalar strain 7 will be 
sufficient. 

For normal crystalline solids there is a rich literature 
spanning decades and centuries of studies of their elastic 
properties. It is well known that for small external strains 
7 normal solids react elastically to create an opposed 
stress a. For very small external strains the stress is 
linear in the strain, with the coefficient of proportionality 
being the shear modulus /x. Choosing an equilibrium 
reference state at 7 = 70 the stress resulting from a small 
strain thus reads 



(T = /i(7 - 7o) 



(1) 



For higher strains the elastic response becomes nonlinear, 
and we denote, for the purpose of the discussion below, 
the nonlinear coefficients as i3„ with n > 2 (that is to say. 



by convention Bi — ii). We thus write for large strains 

f^(7) =m(7-7o) + (7 -70)^ + ^3 (7-70)^ + - •• • (2) 

For perfect monatomic crystalline solids the elastic the- 
ory is entirely transparent since both the linear and the 
nonlinear coefficients in this expansion can be computed 
directly from the change in system energy as a function 
of the strain. 
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(the so-called Born term, and see below for details). 
When the crystalline solid becomes riddled with defects, 
the nature of the elastic theory becomes much more ob- 
scure, since the Born term does not tell the whole story. 
In amorphous solids deviations from the Born approx- 
imation are inevitable. Any strain is bound to bring 
about, in addition to an affine transformation in the posi- 
tion of the particles, also non-affine transformations that 
are necessary to bring back the system to a mechanical 
equilibrium in which the net force fi on every particle is 
zero. Nevertheless, The fundamental assumption that is 
equivalent to the assertion that a piece of material is a 
solid, is that before any plastic deformation takes place, 
the potential can be expanded in powers of the strain to 
derive the various expressions for the stress, shear mod- 
ulus and higher nonlinear elastic coefficients. 

In the following we consider deformations via pa- 
rameterized transformations on the particle coordinates 
J'(7) = I + where h determines the imposed defor- 
mation [3|], 

r,^ J ■r^ + u, , (4) 

where the non-affine coordinates Ui are additional dis- 
placements that guarantee that the mechanical equilib- 
rium constraint is fulfilled. Total derivatives with respect 
to strain in the athermal limit should thus satisfy the 
zero- forces constraint [1, 0, Q : 

d d dui d d dui d 
dj dj d'j dui d"f dj dri 



(5) 
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where the second equaUty results from from Eq. (|4]), and 
here and below repeated subscript indices are summed 
over, unless indicated otherwise. 
Using Eq. ([SJ to compute the stress 
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where the second equality holds because of mechanical 
equilibrium fi = —dU/dri = 0. The conditions of me- 
chanical equilibrium is used further in asserting that the 
amorphous solid remains in mechanical equilibrium also 
after straining (but before any plastic event takes place) . 
Thus 



df^ ^ df^ ^ df^ ^ duj 

c?7 d"f dvj c?7 



0. 



(7) 
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Using the notation for the mismatch forces 
and the non-afline velocities 



(8) FIG. 1: The non-affine velocities v = ^ = -H~^ • H for a 
typical quenched realization of a two dimensional glass form- 
ing model, with iV = 10000. 
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we can rewrite Eq. (O as 



V = -H^ ■ H 



(9) 



(10) 



Here the Hessian matrix H.j^ is defined as iJ,-,- = 



A realization of the non-affine velocities is plotted in 
Fig. m for a two dimensional model amorphous solid. 
With the definition of the non-affine velocities, the total 
derivative operator ([5]) takes the form 
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By repeating the application of the total derivative of 
stress with respect to strain, the shear modulus appears 



1 d^U _ I 
Vdf ~ V 
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where the last step is achieved with the help of Eq. ([T0| . 
In Eq. (|12l) /is is the affine Born contribution to the shear 
modulus while the second term represents the reversible 
non-affine contribution to the shear modulus due to the 
motion of the atoms in the amorphous solid finding new 
equilibrium sites after being strained by an affine trans- 
formation. We should stress that until now the strain 
remains reversible in the sense that reducing 7 back to 
zero (or to 70 if that was the reference state) returns all 
the coordinates to their reference value. 



Examining the expression for the shear modulus, 
Eq. p^ . we note that it contains a contraction of the in- 
verse of the Hessian matrix H~^. The existence of eigen- 
modes associated with vanishing small eigenvalues of the 
Hessian H will result in a singularity in H~^, which in 
turn may lead to singularities in the shear modulus. We 
stress here that we always consider the inverse Hessian 
after removing the Goldstone modes. We will show in 
the following that higher order elastic coefficients consist 
of terms containing an increasing number of contractions 
of H~^, which can thus lead to increasingly higher order 
singularities. 

There are two main physical reasons for the existence 
of low-lying modes in amorphous solids. First, it is well 
known that with increasing system size, one finds low 
lying delocalized modes, for example in the form of plane 
waves with wave- vector of the order of q ~ 1/i, with 
associated eigenvalues A which can be estimated as 
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c q 



-2/d 



(13) 



In addition to the existence of low lying delocalized 
modes, there exist plastic events, which are understood 
as a mechanical instability in amorphous systems 
When a plastic event occurs at some value of the external 
strain 7 = 7p, it is manifested as a saddle node bifurca- 
tion, with the lowest eigenvalue, denoted as Ap, vanish- 
ing like Ap cx y/jp — 7. It was found in [Tj] that when 
N 00 then 7p 0; plastic events occur at smaller 
and smaller values of the strain as the system increases 
in size. This important finding is reiterated below in 
Sect. IIVI This inevitable proximity of mechanical insta- 
bilities is therefore the second reason for the existence of 
low-lying modes, and below we will show that it is in fact 
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the most dangerous one. 

The question we address in this paper is how the shear 
modulus fi and higher order elastic moduli Bn scale in 
the thermodynamic limit N — > oo, V ^ oo and N/V — p 
in athermal amorphous media, by directly analyzing the 
effect of low-lying modes through their corresponding 
statistics. The analytic estimates are then validated 
against numerical simulations of model glass formers, 
both in two and three dimensions. The first point to 
be made is that because of the amorphous nature of the 
quenched solid, the various elastic moduli are random 
variables and are therefore characterized by their expec- 
tation values, their fluctuations, or more fully by their 
probability distribution functions (pdf's). We need to 
know whether these pdf's collapse to a delta function in 
the thermodynamic limit. If not, are the pdf's attaining 
a limit function in the thermodynamic limit, or are the 
fluctuations divergent in this limit. It will turn out that 
/i exists in the thermodynamic limit and its distribution 
attains a delta function. On the other hand, the pdf of 
B2 does not collapse and remains wide, independently of 
the value of N. The coefficient B3 and all higher coeffi- 
cients will be shown to diverge; they do not exist at all in 
the thermodynamic limit. We stress that these consid- 
eration are entirely orthogonal to questions raised, say 
in Ref. P where is was argued that amorphous glassy 
systems at finite temperatures can always flow if we wait 
sufficiently long. The present analysis has nothing to do 
with temperature fluctuations or with long time scales. 
Our findings are also fundamentally different from those 
conveyed in Ref. [1| where it was shown that for a solid 
that contains a crack the series ^ has a zero radius of 
convergence. The present analysis pertains to homoge- 
nous amorphous solid and is a direct quest for the exis- 
tence of its elastic theory in athermal conditions, but in 
the thermodynamic limit. 

The structure of this paper is as follows; in Sect. [Tl] 
we describe the models and numerical methods used in 
this work. In Sect. Illll we derive expressions for the most 
singular terms of the n'th order elastic coefficients. In 
Sect. IIVI we remind the reader when a freshly quenched 
amorphous solid is expected to suffer its first plastic event 
upon increased strain, and introduce the statistics of 
'plastic modes'. Sect. |V] is devoted to the analysis of 
the effect of plane waves on the elastic coefficients in the 
thermodynamic limit; the conclusion is that plane waves 
do not lead to singularities in the elastic coefficients of 
any order. In Sect. |Vl]we present an analysis of the ef- 
fect of localized plastic modes on the elastic coefficients. 
We show that these modes do not lead to the divergence 
of the shear modulus in the thermodynamic limit. How- 
ever, they do lead to enormous fiuctuations in B2 - its pdf 
is predicted to attain a limit distribution that does not 
narrow upon increasing the system size. In addition, we 
predict and exemplify that the plastic modes lead to a di- 
vergence in the mean of B3 in the thermodynamic limit, 
concluding that it does not exist. All the higher order co- 
efficients are expected to diverge even faster. The paper 



is summarized and discussed in Sect. lVlIl We present our 
conclusion regarding the physical meaning of the anal- 
ysis, proposing that our results explain some surprising 
experimental statements about metallic glasses that were 
made by various groups. 



II. MODEL AND NUMERICAL METHODS 

Below we employ a model glass-forming system with 
point particles of two 'sizes' but of equal mass m in two 
and three dimensions (2D and 3D, respectively), inter- 
acting via a pairwise potential of the form 



21 



=0 





(14) 

where is the distance between particle i and j, s is the 
energy scale, and Xc is the dimensionless length for which 
the potential will vanish continuously up to q derivatives. 
The interaction lengthscale Xij between any two particles 
i and j is = 1.0 A, A^ = 1.18 A and A^ = 1.4A for two 
'small' particles, one 'large' and one 'small' particle and 
two 'large' particle respectively. The coefficients 021 are 
given by 



C21 



(-1) 



e+i 



{k + 2q)\\ 



i2q ~ 2£y.\{2£)U (k - 2y.\{k + 2e) 



(15) 



We chose the parameters Xc = 1.48, fc = 10 and q — 3. 
The unit of length A is set to be the interaction length 
scale of two small particles, and e is the unit of en- 
ergy and temperature. Accordingly, the time is mea- 
sured in units of — ^ my? je. The density is set to 
be N jV = 0.85A^^ for our two-dimensional systems, and 
N/V — 0.82A~'^ for our three-dimensional systems. We 
prepared 5000 independent amorphous solids by quench- 
ing high temperature equilibrium states with the rate of 
XO"**^. Any residual heat was removed by a potential 
energy minimization. We deformed our systems using 
the athermal quasi-static (AQS) scheme described in de- 
tail in [Ki] . Elastic coefficients were calculated using the 
prescriptions and microscopic expressions given in 
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MOST SINGULAR TERMS OF ELASTIC 
COEFFICIENTS 



III. 



The shear modulus was calculated exactly in Eq. (|12l) . 
It contains two terms, the Born term and the consequence 
of the non-affine field v. The Born term is always finite, 
being a straightforward partial derivative of the energy 
(this remark pertains to the Born terms yd"'U/d'^" in 
all the higher order elastic coefficients as well). On the 
other hand, the second term in Eq. (IT^ which is due 
to the non-affine response, contains the contraction of 
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which can be potentially singular. The higher or- 
der elastic coefficients have more terms that result from 
the non-afHne response, but it is always easy to identify 
which is the potentially most divergent term. In the fol- 
lowing section, we will consider the expressions for the 
most singular terms of the rt'th order elastic coefficients, 
denoted Z)„, which are the terms containing the largest 
number of contractions with H^. It should be under- 
stood that Bn contains terms of all orders in H~^, from 
the zeroth order (which is the Born term) all the way 
to the 2n — 1 order, cf. 0- It is the latter one which is 
potentially most singular, and we focus on that term in 
the sequel. 

To derive the most singular terms D„ of Bn, one can 
take consecutive total derivatives of the most divergent 
part of /i with respect to the external strain. The techni- 
cal tool to do so is derived by taking the total derivative 
([TT|) with respect to strain of the equation ■ H — I, 
then one obtains the relation ;2] 



dH 



H ' — 



H 



-1 



-H-^-{vT)-H-^ , (16) 



where we have introduced the third-rank tensor T^fe = 



drkdr Or " Sincc V contains one contraction of , cf. 
Eq. (fTO|) , the second term on the RHS of Eq. ([T6l) contains 
a total of three contractions of H^^, as opposed to only 
one contraction in the first term. So it is the second term 
■ {v ■ T) ■ that gives rise to the most singular 
terms which are considered below. 

Using Eq. (|16p. we can operate on the most singular 
term of the shear modulus, 



V 



(17) 



to obtain the most singular terms Z)„ of the elastic coef- 
ficients Bn of the next two orders: 



Do = 



Twvv 



V 



and 



Xvv : T) ■ ■ (T 

V 



(18) 



(19) 



The full tensorial form of D2 and was worked out in 
, and one can check directly there that Eq. ([T9| indeed 
presents their most singular terms. 

We denote the set of Nd eigenvalues of the Hessian 
by {^k}k=i and their associated real eigenvectors by 
|\I>(*:)|^^ For concreteness, we denote delocalized plane 
wave modes as and their associated eigenvalues as A. 
Similarly we will denote localized plastic modes as 4*, 
and their associated eigenvalues as A. We next introduce 
the contractions and bkim- 



a-k 



J' '■ \f (fe)vi>W\i>('^ 



(20) 



Using the normal mode decomposition of the inverse Hes- 
sian, = J2k ^' x^'^ ' 1 we write the most singular 
terms employing the contractions and bkim as 



and 



3 
V 



Do 
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ijkim 



-E 



aka£a,nbkem 
AfcA^A m 



k£m 



(21) 



(22) 



The n'th most singular term can be written in a symbolic 
notation as 



n - 



Cn \ - {a} 
V ^ 



n+1 



{6}"- 



(...) 



{A}2»- 



(23) 



with the combinatorial factors c„ = (2n — 3)!!, and the 
sum should be understood as running over the corre- 
sponding 2n — 1 indices of the eigenvectors. We list a 
number of important realizations regarding Eq. ([23^ : 



(z) Every index that is being summed upon appears 
exactly twice in the numerator with a correspond- 
ing eigenvalue in the denominator. This is a direct 
result of the contraction of the normal mode de- 

-1 ^('^)^('^) 



composition of H 



Ek 



(a) The sum can consist of only one arrangement of 
indices (up to trivial permutations in identities of 
dummy indices). This can be understood by ex- 
amining the form of Eq. (|16|) : in the various con- 
tractions, any is always contracted with two 
tensors, never with a single one. So, for instance, 
contractions of the form ^ : H~f} cannot appear. 
Consequently, combinations of the form bikk or bkkk 
cannot appear. 



(in) For odd n, I?„ is can be written as 



Dn = - 



where 



Sfc = E 



(...) 



{A}- 



(24) 



(25) 



and the above sum should be understood to run 
over n — 1 corresponding indices. This is a con- 
sequence of the form w ■ ■ w that the most 
singular terms of odd order take, as seen in, e.g., 
Eqs. (dZl) and (HSl). 

The above points are consistent with the symmetry to 
flipping 7 — > —7 that require that for even n (Dn) = 0. 
One should note that although i?„ are defined as the n'th 
total derivatives of stress with respect to strain, the most 
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singular terms contain at most a third derivative of the 
potential energy with respect to coordinates. This stems 
from the fact that the most singular terms contain the 
highest number of contractions with H~^\ dimensional 
consideration then limit the order of the energy deriva- 
tives with respect to coordinates to third order. This is 
the order that contributes since the first order derivative 
vanishes due to mechanical equilibrium, and the second 
order is the Hessian itself which annuls . Of course 
one finds also mixed second order derivatives with respect 
to external strain and coordinates appearing through H. 



IV. DISTRIBUTION OF PLASTIC MODES 

To help us in studying the issues raised in this pa- 
per, we briefly review one of the findings of Ref. Q, that 
plays a key role in the analysis of the thermodynamic 
limit throughout this work. Ref. [3] employed the same 
generic glass-forming system in two and three dimensions 
as introduced above, see [7] for details. The numerical ex- 
periments performed were as follows: an ensemble of in- 
dependent initial glassy configurations quenched from the 
high temperature liquid was constructed. Then the AQS 
scheme (see [ICi] for details) was utilized to strain each 
system up to the first mechanical instability occurring at 
some strain value A7iso, see Fig. [2] for an illustrating car- 
toon. The mechanical instability is a saddle node bifur- 
cation in which the lowest non-zero eigenvalue Ap of the 
Hessian matrix vanishes like Ap cx VTp^^- Statistics 
of A7iso were collected for a variety of system sizes. The 
first plastic event (when a freshly quenched un-deformed 
system is strained) does not occur for any infinitesimal 
value of 7 and careful measurement of the mean strain in- 
terval (A7iso) that separates the un-deformed state from 
the first plastic event results in a scaling law 

(A7i,o) ^ , Aso ~ -0.62 . (26) 

In Fig. 13] we show the numerical evidence that supports 
Eq. If we accept (the reasonable but not proven) 

assertion that the scaling law measured at relatively low 
values of N continues for all iV, this scaling law means 
that the potentially unstable eigenvalue of the Hessian, 
Ap, which goes through the saddle- node bifurcation at 
the first plastic event scales with the system size like 

Ap^ v/A^^iV^'-/^ (27) 

We should state at this point that the accuracy of mea- 
surement cannot rule out that the second digit in the 
exponent is not accurate, so we need to use Eq. (|27|) 
with a grain of salt, as discussed below. Based on this 
finding, we work out the accumulation of localized modes 
around Ap. Once a distribution of 'plastic modes', -P(A), 
is formed in the thermodynamic limit, we only know that 
when we observe its minimal value we find Ap - 7^^1-/2. 
Using the WeibuU theorem fT2!| we can therefore estimate 
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7 7p 



FIG. 2: Cartoon of a typical stress vs. strain curve of a single 
instance of the numerical experiment performed in [J. Each 
un-deformed system was strained until the first mechanical 
instability was encountered at some strain value 7p. Statistics 
of A7iso = 7p were collected for a variety of system sizes, both 
in two and three dimensions. 
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FIG. 3: The mean strain interval before the first mechanical 
instability (A7iso), for 2D (left panel) and 3D (right panel). 
The continuous lines represent the scaling law Eq. (|26p . 



that 

P(A) ^ {N/\d)\^K{\) , (28) 

where 6 — — 2+/^'°° ^ 2.23, and K{\) is some smooth 
function with a cutoff of the order of a typical (Debye) 
cutoff eigenvalue \d- 



V. EFFECT OF DELOCALIZED MODES ON 
ELASTIC COEFFICIENTS 

In this section we show that the delocalized modes in 
the form of plane waves do not contribute to the diver- 
gence of the elastic coefficients of any order. A reader 
who is more interested in the divergence due to the plas- 
tic modes can skip this chapter in favor of the next one. 

Before we turn to the analysis, we work out the system 
size dependence of the contractions Eq. (|20p when calcu- 
lated in the context of delocalized plane- wave modes 
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Scaling of the contractions a and b on 
delocalized modes 



with Wij = Wj — Wi. Using this relation and Eq. (|30l) . 
we can estimate 



We begin with cife = H • 

df 



97 



where fij is the force exerted by the jth particle on the 
ith particle. Using the symmetry of the binary force we 
can rewrite 



2 ^ 



(fc) ^ 



(fe) 



9 ' 



(29) 



where we have introduced the notation '^1''^= 'i'l'^ — . 
For long wavelength modes 'S'^'^'' ~ sin(qfc • r^) and 
we expect that the maximal contribution to the difference 

~ (k) 

will follow the scaling 



(30) 



Since can be positive or negative, the sum of N 

terms in Eq. (|29p cancels the y/N of the normalization of 
the eigenmode and accordingly 




N ^ —\/XkXi\„ 



(34) 



where the factor of vTV accounts for summing 0{N) ran- 
dom terms which can be positive or negative. 



B. Proof of no divergence 

With the scaling laws pip and ([M)) for the contractions 
a and h respectively in hand, we are in the position to 
calculate the contribution of delocalized modes on the 
elastic coefficients. 

We begin with the case of odd n; we return to the 
scaling of Sk^ given in a symbolic notation in Eq. (j25p . and 
consider the sums as running over the delocalized modes 

Since {a}T^{6}^ - \/a^{A}"- VTV"^ , then Sk is 
a sum over ©(TV"^^) positive or negative terms of order 

V^{A}"~V^^- This leaves us with Sfc \/Yk. Next, 
plugging in this scaling into p4)). we obtain 



y^Xk 



Oil) , odd 



(35) 



since it is (the negative of) a sum of N positive terms, 
each of 0(1), and the pre-f actor oi 1/V cancels the N 
dependence altogether. To exemplify these arguments 
for odd n, consider D3 given in Eq. (j22p : we can rewrite 
it as 



a-k ^ <lk 



(31) 



We turn now to the contraction of delocalized 
plane-wave modes on the third-rank tensor bkim = 

Tm^:*^'"^**-^^*^™''- To proceed we need to use the sym- 
metry properties of the tensor that were spelled out in 
Including spatial component Greek indices in the 
superscripts, the symmetries are 



rpyrix 
^ijk 

ijj 

Hi 



iij 

^ijk 



= 



iij 

ErpVVX 
iij 

rpl^Xn _ rpXt^V 

iji jii 



if i^j^k 

if i^j 



if i^j 



(32) 



'^ijk' ^^'^ poss. perm, of vrjx 



Using these symmetries, it is easily verified that for gen- 
eral vectors w,y,z: 



3 1 didjbijk agambkim 
^ J, Afe ■■ XiXi XfX„ 



'^''1 em. 



3 \ ^ 1 j \ ^ aidjbijk 



In this case didjbijk ^ XiXjy Xk/N and 



Sk 



didjbijk 

ij 



XiXj 



(36) 



since it is a sum of 0{N'^) positive or negative terms, 
each of order \/Xk/N . From here, the scaling given in 
Eq. (p5)) is immediately obtained. 

We next consider the case of even n; starting from 
Eq. (j23l) . we again consider sums over the delocalized 
modes ^ . In the even n case we consider the numerator 
~ {A}2"-77V"-i, then 



^Tijk-w,yjZk = 



ijk 



2 Tiij.WijUijZij 

i,j=^i 



(33) 



{A}2"-i 



1 



(37) 
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Now Dn is a sum over 0{N'^"^~^) positive or negative 
terms of order 1/iV" (incorporating the 1/V pre-factor), 
hence, for even n 



1 



N 



even n 



(38) 



To exempHfy these arguments for even n, consider Dn 
which can be worked out from Eq. 



Da 



15 

17 2^ 



(39) 



In this case, using again ([51]) and (IMl) . we have 

with a sign which can be positive or negative. Therefore, 
the sum in (p9| runs over 0{N'^) terms, each of order 
1/7V^. This, together with the 1/V prefactor, results in 
the scaling D^ ^ 1/ Vn, as expected from Eq. (pS)) . 

We conclude that the effect of low-lying delocalized 
modes on the most singular terms of the elastic coeffi- 
cients is regular; we obtained the scaling Z?„ ^ £>(!) for 
odd n, and D„ ~ for even n, which is consistent 

with the symmetries and the intensive nature of the elas- 
tic coefficients. 



VI. EFFECT OF PLASTIC MODES ON 
ELASTIC COEFFICIENTS 

The most important physical origin for low lying 
modes, besides the presence of delocalized modes, is the 
eminent proximity of plastic failures manifested as me- 
chanical instabilities. In this section we analyze the effect 
of the resulting localized plastic modes on the elastic co- 
efficients of various orders. As opposed to the case of 
delocalized modes (c.f. Sect. [V]), in this case we cannot 
present an a-priori estimate of the system size scaling of 
the contractions Eq. ([20l) when calculated for the case 
of plastic modes 'S'. We will need to rely on numerical 
information. Notwithstanding, we go as far as possible 
with scaling arguments before turning to the numerical 
simulations. 

We define the contractions for the case of localized 

plastic modes, bki^ = T : 'i'C^)*^*^'") and = H • 
'S'^'^'. Since we expect the plastic modes to be highly 
localized, the triple contractions bkim can be taken to 
be represented approximately by delta-functions bkem ^ 
SkiSim- With this approximation (whose consequence we 
test below against our numerics) we conclude that all the 
sums appearing in Eq. (|23)) (when considered to run over 
the plastic modes '4'), are dominated by the diagonal 
terms. 
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FIG. 4: An example of a of a single measurement of D2 close 
to a mechanical instability at 7p in 3D for A'^ = 1000. We 
extract the coefficient of the resulting power law D2 ~ (7p — 
7)" 2 J and average over realizations to obtain Fig. [5] 



It has been shown 0, that close to a mechanical 
instability, the eigenvalue Xp associated with the eigen- 
mode which is becoming unstable follows the equation of 
motion 



dXr 



c?7 



-apbpX, 



O(Api) 



Sfefefe 



for k 



: P . 

(41) 

This equation is analogous to Eq. ([T5|. where it is given 
in terms of . The analogous equation for A^^ is more 
involved, but once we apply the above assumption that 
bkem ~ SkeSem it simplifies to an isomorphic equation 



dXl 



c?7 



-akhXf + OCK') 



(42) 



A consequence of this result is that the product aubk has 
to remain independent of the system size. We see this 
from Eq. which implies that the LHS of Eq. (gl]) is 
of the order of X^^ . 

The assumption that akbk is system size indepen- 
dent does not guarantee that a^~^^b'^~^ is also system 
size independent. In order to establish this, we can 
measure, say, D2 close to a mechanical instability at 
7p. Defining ap and bp as dk and b^ for the plastic 
mode '4'(^\ such a measurement results in a power law 
D2 — (apSp/y)(7p — 7)" I [2I, allowing us to extract 
the pre-factor. By taking statistics over the system size 
dependence of the pre-factor, one can then deduce how 
d^^^b^~^ scales with system size. 

We now turn to the numerical simulations. For every 
independent system from our ensemble of 5000 realiza- 
tion per system size, we strain the system using the AQS 
scheme, until the first mechanical instability is encoun- 
tered at 7p. We then backtrack and measure the most 
singular term D2 at various distances from the instabil- 
ity 7p — 7- An example of a single such measurement is 
plotted in Fig. |4l We then extract the pre-factor of the 
resulting power law, and collect statistics of this number 
for various system sizes, both in two and three dimen- 
sions. The results are plotted in Fig. [S] We find that the 



8 



10' 



a, 

<C3 



10 



(2D) 



□ □ □ 



10' 



10 



(3D) 



o o o o o o 



10 



N 



10 



10 



N 



FIG. 5: The mean of the product of the contractions a'^bp, 
as a function of system size A''. For the calculation below the 
weak observed dependence is taken as no dependence ap ~ 
N°. 



dependence of (ap6p) on N is very weak. Together with 
the assumption that apbp is N independent we conclude 
that to a reasonable approximation we can take 



(43) 



both in two and three dimensions. We will assume in 
the following that all the a| have this N independence, 
without a systematic dependence on the index fc, i.e. 

Having established the scaling of a^, we next consider 
the sum in Eq. (1401) : we begin with n = 1, and write this 
sum as an integral over the density of plastic modes P{X) 
given in Eq. then 
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N 



X'^K{X)dX 
X 



(44) 



N TT^ 



Since sa 2.2 > 1, it is obvious that the low- lying plastic 
modes do not lead to a singularity in the most singular 
term of the first order modulus Di. 

For an odd n > 1 the sum in Eq. (|40l) becomes insen- 
sitive to the alternating signs of akbk', 
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n+1 Tn— 1 
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\2n-l 
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From here 



E 



2n-l 



X^K(X)dX 
A2"-i 



(45) 



(46) 



N TT? 



For an even n the sum in Eq. (|40l) consists of positive 
and negative terms, and whether the sum diverges or 
converges depends delicately on the exponents Q and 2n— 
1. Since the mean of the moduli with even n vanishes, we 
are interested in their pdf 's (which are always symmetric 
around zero). The issue of convergence or divergence is 
dealt in detail in Appendix [XJ where we study the pdf 's 
of the sum 



(fe' 

.9± 



(47) 



where g\. is -1-1 or -1 with equal probabilities. We show 
that for the value Q « 2.2 and n > 1 all these pdf's tend 
in the thermodynamic limit to a form with power-law 
tails. The result is that the re-scaled pdf 



/ 



2 71-1 



/(y) ~ \y\ 



1 , n even, \y\ > 1, 

(48) 

2T1-1 

meaning that once the pdf's are rescaled by N i+» they 
attain a limit form. 

We can now finally turn to estimate the system size 
scaling of D„; incorporating the results of Appendix 1X1 
with Eqs. (|ini),(|lSI) and dH]), we estimate 



1 2,.-l 
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For n = 2, the prediction is 



(49) 



(50) 



Note that with our numerical accuracy regarding the ex- 
ponents we cannot rule out that the exponent is actually 
zero. We expect therefore the distributions of Z?2 to be 
either independent of system size or to decay very slowly 
to a delta- function. For n = 3, we calculate 



0.56 



(51) 



Consequently (-D3) is predicted to diverge in the ther- 
modynamic limit, leading to the divergence of the third 
order elastic coefficients B3. Needless to say, higher or- 
der coefficients will display stronger divergences, and will 
not exist in the thermodynamic limit. 

Turning again to the numerics, we measured the most 
singular terms Di,D2 and D3 for each member of our 
ensemble of quenched amorphous solids, both in 2D and 
3D. For each member of the ensemble we computed these 
objects directly from their definitions ([TT]) - (IT51) . Due 
to the amorphous nature of our solid, each member of 
the ensemble provides a different value to the measured 
quantity. The distributions of Di are displayed in Fig. |51 
as predicted, they show no anomalies, and their width 
decay with increasing system size, as expected from an 
intensive thermodynamic variable. 

Next, we present the distributions of D2 in Fig. [T] As 
required from isotropy, the distributions are symmetric 
around zero, and (-D2) — 0. As predicted, the width 
of the distribution remains N independent (or decreases 
extremely slowly) , meaning that the fluctuations over re- 
alizations of D2 do not decay appreciably with increasing 
the system size, in spite the intensive nature of the elastic 
coefficients. 

In Fig. [5] we display the distributions of D3 measured 
in our simulations. Indeed, we find that (-D3) grows with 
system size, in agreement with the estimate given here. 
In Fig. |9]we show the distributions of D3 rescaled by the 
expected A^-dependence given by Eq. the data 

collapse is in very good support of the estimated expo- 
nents. We measured approximately (^20 ~ 0.62 in 2D 
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FIG. 6: Distributions of Di measured for each instance in 
our ensemble of quenched amorphous sohds. As expected, 
the width of these distributions decays with increasing system 
size. 
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FIG. 7: Distributions of D2 measured for each instance in our 
ensemble of quenched amorphous solids. As expected, the 
width of these distributions does not decay with increasing 
system size. 



and CsD ~ 0.56 in 3D. Note that all the elastic coeffi- 
cients from i?3 and higher are expected to diverge in the 
thermodynamic limit even faster. 
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FIG. 8: Distributions of D3 measured for each instance in our 
ensemble of quenched amorphous solids. 
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FIG. 9: Distributions of Dz rescaled by the estimated N- 
dependence, cf. Eq. (|5ip . The data collapse is in excellent 
agreement with the predicted scaling behavior. We find the 
best collapse with C,2d ~ 0.62 for the 2D data and C,zd ~ 0.56 
for the 3D data. 



VII. SUMMARY AND DISCUSSION 

We have learned in this paper that in amorphous solids 
of the type discussed above the elastic theory exists only 
for infinitesimal strains in the thermodynamic limit. We 
have found that the shear modulus exists in the ther- 
modynamic limits, but starting with B2 the pdf's of all 
the nonlinear coefficients i?„ do not converge to a delta- 
function in the thermodynamic limit. Starting with 
all the coefiicient diverge with N . The final result can be 
summarized by the prediction 



(52) 



Note that the exponent 9 is not expected to be universal 
since it depends on the quenching protocol and other 
details of the system. We expect however quite generally 
that 9 > 2 [1^. Thus the shear modulus exists always, 
providing us with a small linear regime where elasticity 
theory may remain useful. 

We reiterate here that the analysis presented here has 
nothing to do with time-scales. As long as our straining 
rate is slow enough, 7 <C i/c the formation of an affine 
transformation is co-temporal with the creation of forces 
on the particles that result in a concurrent non affine 
transformation which is the source of potential singulari- 
ties in the theory. One cannot 'freeze' the non affine part 
of the transformation and go on straining. The 'velocity' 
V above is real, and represents how the non affine trans- 
formation occurs on the same time scale that the strain 
is increased. 

The upshot of the analysis presented above is that in 
amorphous solids in the thermodynamic limit it is im- 
possible to separate elastic from plastic responses. They 
are mixed together in an intimate way, defying the exis- 
tence of a naive nonlinear elastic theory in the form of a 
Taylor expansion. It is possible that one needs to renor- 
malize the theory in order to obtain a well posed expan- 
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sion in another variable; this possibihty will be explored 
in a forthcoming publication. At this point our feeling 
is that the results presented above explain the surpris- 
ingly high degree of plasticity that was found in recent 
experiments on metallic glasses [13; ] , and of the fact that 
macroscopic samples of metallic glasses cannot be main- 
tained in the elasto-plastic steady state that is found in 
small sample simulations [14 , 15] : they appear to break 
as soon as they reach the yielding transition. Also, we 
point out to the interesting experiments [H, [l3| where 
it was shown that in nano-samples of metallic particles 
one could reach the elasto-plastic steady state. While 
not mathematically related to the analysis shown in this 
paper, we propose that all these observations are in ac- 
cordance with the theoretical conclusions presented here. 
It may be worthwhile to examine systematically the sys- 
tem size dependence of the existence of the yielding tran- 
sition to a steady elasto-plastic state in metallic glasses 
and other systems of the type presented here. In par- 
ticular, it should be interesting to examine the nonlinear 
elastic coefficients in small system, like nano-samples of 
metallic glasses, and attempt to confirm the divergence 
of i?3 and of higher order nonlinear elastic coefficients. 
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Following this formulation, we ran simulations of the 
above statistics and extracted the pdf's of for n = 2 
and n — A, which are the cases of interest in Sect. I VII 
The pdfs for X2 and X4 are plotted in semi-log scales 
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FIG. 10: Distributions of X2/NTTe^ calculated using 9 
see text for details of statistical model. 



2.2, 



in Figs. [TU] and [TTJ As expected from the symmetry of 

(k) 

the factors g^. , the distributions are symmetric around 
zero. We confirm empirically the known result [l9j that 
the pdfs of Xn{N) for various iV's when plotted as a func- 

2n-l 

tion of the re-scaled argument Xn/N i+f , attain a limit 
distribution, known as a 'Levy law' [l9| . Our simulations 
were carried out using 9 = 2.2, see Sect. IIVI and Sect. IVll 
for the origin of this exponent; we note that the extracted 
scaling law does not apply to any choice of exponents 9 
and 2n—l, see review by Bouchaud and Georges [l^ for 
further details. According to [l^, one expects the lim- 
iting distributions the sums X„ for higher orders n > 4 
to obey the same scaling, namely they are given by a 

2ra-l 

scaling function of Xn/N 1+" . 



Appendix A: System size scaling of ^^(A'^) 

The sum Xn{N) was defined in Eq. (|T7|) : here we pro- 
vide a numerical analysis of the statistical problem, from 
which we will deduce the scaling of Xn with the system 
size N . For clarity, we repeat the formulation of the 
problem: given a probability distribution K{X) ^ for 
small A, we draw 0{N) numbers from this distribution, 
and calculate 
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where as before g[^^ = ±1 with probability 1/2 for each 
possibility. Notice that here we only consider the case of 

even integers n, then (.9±^'' ) = g'±^ ■ 



FIG. 11: Distributions of Xi/N^^ , calculated using 
see text for details of statistical model. 



6 = 2.2, 
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